On the Wake Structure in Streaming Complex Plasmas 



Patrick Ludwig 1 , Wojciech J. Miloch 2 ' 3 , Hanno Kahlert 1 , and Michael 
Bonitz 1 

Institut fur Theoretische Physik und Astrophysik, Christian-Albrechts-Universitat zu Kiel, 
24098 Kiel, Germany 

2 Department of Physics, University of Oslo, Box 1048 Blindern, 0316 Oslo, Norway 

3 Department of Physics and Technology, University of Troms0, 9037 Troms0, Norway 

E-mail: ludwigOtheo-physik .uni-kiel . de 

Abstract. The theoretical description of complex (dusty) plasmas requires multiscale concepts 
that adequately incorporate the correlated interplay of streaming electrons and ions, neutrals, 
and dust grains. Knowing the effective dust-dust interaction, the multiscale problem can be 
effectively reduced to a one-component plasma model of the dust subsystem. The goal of the 
present publication is a systematic evaluation of the electrostatic potential distribution around 
a dust grain in the presence of a streaming plasma environment by means of two complementary 
approaches: (i) a high precision computation of the dynamically screened Coulomb potential 
from the dynamic dielectric function, and (ii) full 3D particle-in-cell simulations, which self- 
consistently include dynamical grain charging and non-linear effects. The applicability of these 
two approaches is addressed. 



1. Introduction 

Complex (dusty) plasmas have been proven to be an instructive reference system for strong 
coupling and correlation effects [IJ |2] |3] . Thereby, complex plasma research inherently connects 
key issues from several fields including low temperature physics, warm dense matter physics, 
surface and solid-state physics, as well as material science [H [5j [U [7] [EJ [9] . 

The first principle description of complex plasmas is a theoretically challenging problem. A 
dusty plasma consists of electrons, positive ions, neutral atoms, and micrometer-sized dust grains, 
i.e., components with distinct mass asymmetry giving rise to dynamics on very different spatial 
and temporal scales. This multi-scale behavior, combined with the non-ideal character of a 
complex plasma, makes a full temporal resolution on all scales a numerically unfeasible task with 
current computer technology. However, there are two numerical approaches that can address 
the challenges associated with the complex plasma system: the One-Component Plasma (OCP) 
model and the Particle-In-Cell (PIC) method. 

The OCP approach relies on the idea to reduce the quasi-neutral, multi-component complex 
plasma to the subsystem of dust grains which interact via a screened Coulomb potential. 
In particular, the statically screened (Yukawa-type) potential yields good agreement with 
the experiments for various specific setups, see for example Refs. [10, llj. Also, this 
isotropic potential is well suited for fundamental analytical and numerical investigations on 
cooperative effects such as self-organized structure formation j3j [12], collective dynamical 
processes |13| [TH II 5 j . spectral properties |16| 117] . or the melting behavior |18l [19] in strongly 
correlated (screened) Coulomb systems. The multi-component and non-equilibrium nature of a 



complex plasma requires, however, a careful analysis of plasma streaming and dynamical grain 
charging effects, which were experimentally shown to become of utmost importance at least 
near the plasma sheath, see e.g. [20 ] |2T| [22] [23] [24] 125] . In a plasma with ions streaming at 
a uniform velocity, the dust potential becomes (strongly) anisotropic and takes the form of an 
oscillatory wake structure in the downstream direction, which was subject to various analytical 
and numerical studies e.g. [26l Ell [28l [29l [30l [Ml E21 [Ml [Ml This 
wake can result in attractive, non-reciprocal forces between the equally charged dust grains and 
lead to remarkable structural and dynamical consequences |43^ 0H SSJ 06J H7J 0U S3 ED 152] . 
Streaming effects can be incorporated in the OCP dynamical screening model by means of Linear 
Response (LR) theory, as presented in section [2] 

The second powerful numerical ansatz to tackle the multi-scale problem of the different plasma 
components is the PIC approach, in which the trajectories of plasma particles are followed in 
electric fields self-consistently [55J. To obtain a self-consistent charge on the grains, the plasma 
dynamics need to be resolved in the vicinity of the grain on spatial and temporal scales which 
are smaller than the Debye length and electron plasma period, respectively. However, due to the 
large mass of grains (as compared to electrons and ions) as well as the large distances between 
the grains, the study of the dynamics of several dust grains requires simulation times that are 
much longer than the dust charging (which is on the ion plasma period time scale). Thus, a self- 
consistent 3D PIC calculation, i.e., a dynamic evolution of the 3D system with more than one 
dust grain with a temporal resolution on the electron or ion plasma period time scale, required 
to ensure self-consistency, easily exceeds the capabilities of present high-performance computing 
systems (5TJ [56| . 

Consequently, short-time, small-scale PIC simulations of dust charging or, alternatively, less 
demanding LR calculations of the dynamical plasma shielding have to be coupled with large-scale 
OCP simulations which incorporate the interaction between the grains on a more abstract level 
|46l H71 148j . Such a coupled multi-scale numerical approach may ensure the description of the 
correlated system dynamics with proper charges on the grains as well as an accurate potential 
distribution in the vicinity of grains. With the goal to form a method-spanning picture, the 
present publication comprises a systematic analysis of the electrostatic potential distribution 
around a dust grain in the presence of a streaming plasma environment by 

(i) a high precision computation of the dynamically screened Coulomb potential from the 
dynamic dielectric function, and 

(ii) a critical assessment of these LR results, in particular in the view of non-linear effects and 
dynamical grain charging processes by means of self-consistent 3D PIC simulations. 

The discussion of the LR results, section [2] extends over a broad range of plasma streaming 
velocities, different electron-to-ion temperature ratios and includes an evaluation of the influence 
of collisional damping. The corresponding 3D PIC simulations, section [3j are performed for 
the collisionless limit and apply, therefore, to the low pressure limit. Finally, it is vital to 
verify whether and to what degree the results obtained by a linear superposition of single grain 
potentials do coincide with the corresponding full 3D PIC simulations. Such a direct comparison 
is performed for the case of two grains in section [4] 

1.1. Plasma Parameters 

The plasma parameters are scaled to relevant base units, which are, in our case, the Debye 
length Xo a = y/sok^Pa] {n a q^) and the plasma period r a = 2Tr/u) Pa , with the corresponding 
plasma frequency uj pa = y/n a q^/(eQm a ), where a denotes electrons (e) or ions (i). T a is the 

1 For a review of previous experimental and theoretical work on the structure of the dust (wake) potential, we 
refer the following in-depth references [5, 6 7 , 23 53 54]. A recommended critcal discussion of previous theoretical 
work on the structure of plasma wakes can be found in |42| . 



corresponding temperature, n a density, m a mass, q a charge, kj$ the Boltzmann constant, and 
Eq is the permittivity of vacuum. The relative streaming velocity Uj of the ions with respect to 
the dust is characterized by the Mach number M = Ui/c s , where the ion sound (Bohm) speed is 
given by c s = ^Jk^T e jmi. Furthermore, the thermal velocity is VT a = \/kBT a /m a . 

The dynamically screened dust potential and in particular the characteristics of the wake, are 
studied with respect to three dimensionless plasma parameters: 

• the ion drift velocity, expressed in terms of the Mach number M, 

• the ratio of electron-to-ion temperature T r = T e /Ti, 

• the ion-neutral scattering frequency (normalized to the ion plasma frequency) Vi n juj Vi . 

While the temperature ratio T e /T, controls the effect of Landau damping, the frequency 
ratio Vi n /uj Vi determines the effect of collisional damping and can be effectively controlled in 
experiments by changing the neutral gas pressure. A simple approximation for the pressure as 
function of the ion-neutral collision frequency and ion temperature is given by 

_ Vin/upj ■ eo^k B Tj ■ n e /e 

0~in 

where eo = 1.602 • 10~ 19 C is the elementary charge, and <7 m ~ 5.0 • 10~ 19 m 2 is the ion-neutral 
collision cross-section as defined in Ref. [57] J^] 

For the LR calculations, we consider Argon, which is a typical working gas in complex plasma 
experiments |5]. The ions are considered to be singly charged (gj = eo). Consequently, ion 
and electron densities are equal, and are set to n e = rii = 2.0 • 10 14 m~ 3 and are considered 
as spatially homogeneous distributions. The considered electron temperature is T e = 2.585 eV 
(~ 30000 K). The corresponding electron Debye length is \[) e = 845 /im. An atomic mass 
of Argon m n = mi = 6.634 • 10~ 26 kg leads to the Bohm speed of c s = 2500 m/s. For the 
considered set of parameters, the ion plasma frequency takes the value of u Pi = 3.0 • 10 6 Hz. For 
the PIC simulations, we simulate ions with a reduced mass mi/m e = 120, and thus uj pi and c s 



are increased accordingly (see section 3.1 for details). The LR and PIC results presented in this 



paper do not involve any kind of free fit parameters. 



1.2. Grain charging and wake formation 

The dust grains are charged by plasma currents and other currents, such as induced by the 
photoelectric effect or secondary electron emission. For the charge equilibrium, the net current 
to the grain is zero, and such a grain is at floating potential with respect to the surrounding 
plasma. If charging is only due to plasma currents (i.e., electron and ion fluxes to the surface), 
then in electropositive plasmas, the charge Qd on the grain and the floating potential will be 
negative, due to the high mobility of electrons. The floating potential on a spherical, finite-sized 
grain can be approximated with the Orbit-Motion-Limited (OML) approach for both stationary 
and flowing collisionless plasmas |51j . We note that in the presence of collisions, the OML theory 
can give incorrect results, as it has been noted by several authors |59[ l60"| 161) . The presence of 
trapped ions created by ion-neutral collisions can increase the ion current to the grain, and thus 
reduce the charge on the grain |62l 163} [64"] . 

In a stationary plasma, the charging of a grain takes approximately one ion plasma period, 
n, a time in which the ions can collectively respond to the perturbations in the electric field. A 
typical charging characteristic for a single, spherical conducting grain is shown in Fig. [TJ where 
the data is taken from 3D PIC simulations for the stationary plasma with T e /Tj = 30. In the 

2 We should note that the ion-neutral scattering frequency Ui n and not the roughly approximated pressure p is 
used as input parameter for the LR calculations. For the energy dependence of the ion-neutral cross-section in 
Argon see e.g. |58j . 
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Figure 1. Typical charging characteristics for a spherical conducting grain in stationary 
plasma. Data is from 3D PIC simulations of a stationary plasma with T e /Tj = 30, time t is 
normalized with the ion plasma period Tj. 



static case, the charge on a small grain can be calculated with the capacitance model for a given 
floating potential [65j: 

Q d = 47re a (l + ^ J $ fl , (2) 

where a is the grain radius and \d the total Debye length defined in terms of electron and ion 
Debye lengths A^ 2 = A^ 2 + A^, 2 for a static case, i.e., M = 0. When charging is only due 
to plasma currents, the floating potential in a collisionless plasma is typically on the order of 
~ —2k^T e /eo. Thus, from the capacitance model one can derive a formula for an approximate 
charge on a negligible small grain with a <C A/?: 



Qd/e = -U00a [lim] T ej[eV] , (3) 

where the grain radius a is given in micrometers, and temperature in eV. In the case of no plasma 
flow relative to the dust grain, the negatively charged grain is shielded by a stationary, positively 
charged ion cloud. The electric potential around a small (non-absorbing) charged grain takes the 
form of a statically screened Coulomb potential which is typically referred to as Debye-Hiickel 
or Yukawa potential 

*YukW = j^-— e~& . (4) 
Aireo r 

Plasma absorption, ion-neutral collisions, and plasma flow can modify the potential distribution. 
In particular, the long-range asymptote of the plasma-mediated potential is often not screened 
exponentially, but has a power-law decay instead, see e.g. |33l 153"! 154") Wf\ 168] . 

The plasma flow breaks the symmetry in the grain charging. In the mezothermal velocity 
regime, when the flow velocity u\ is larger than the ion thermal velocity and lower than the 
electron thermal velocity vt £ (^T; < u i < v T e ), the grain will receive a larger ion flux to the 
surface on the upstream than on the downstream side, while the electron flux will be similar to 
both sides. Anisotropic charging can lead to an electric dipole on the grain [30, 66]. The electric 



dipole moment depends on the material of the grain and is more pronounced for supersonic 
flows [37J. For small conducting grains the electric dipole can often be neglected, as the charge 
is redistributed on the surface, but for larger conducting grains, potential distribution in the 
vicinity of a grain can modify the charge distribution on the grain surface. 

A streaming plasma leads to a wake in the density and potential distributions around the grain. 
The main mechanisms for the wake formation are plasma absorption on the grain surface and 
the influence of the electric field on plasma particle trajectories in the vicinity of the grain |41] , 
Plasma absorption leads to density depletion on the downstream side, while electric fields can 
scatter ions into the wake region forming the ion focus region. Both effects are more pronounced 
for cold ions, as well as for supersonic flows, when also a Mach cone forms. 

The charging time in a streaming plasma is longer than in a stationary plasma and depends 
on the material of a dust grain, being longer for insulating grains, for which it can take several 
ion plasma periods to reach stationary conditions. 



2. Linear Response Approach 

According to the OCP concept, the multi-component dusty plasma can be reduced to the 
subsystem of dust grains which interact via a dynamically screened Coulomb potential. This 
dynamical potential takes into account the effect of plasma streaming on the dielectric response 
of the plasma. Depending on the quality of dielectric function e(k, oj), the considered dynamical 
screening model comprises an accurate representation of the most essential plasma properties 
including screening, wakefield oscillations, ion and electron thermal effects, Landau damping, as 
well as collisional damping |33| 08] . 



2.1. Theoretical Approach 

Considering a weak (linear) response of the plasma to the presence of the dust grain, the electric 
potential can be computed by a 3D Fourier transformation 

Qd 
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where r denotes the distance from the grain's center of mass. The (shifted) Maxwellian 
plasma is represented by the following longitudinal dielectric function which includes 
Bhatnagar-Gross-Krook (BGK)-type ion-neutral collisions |69| 170] 
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where we use the substitution 
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The screening by electrons is considered as static (Yukawa-type), since the electron flow is of no 
relevance (u e <C vt e )- The plasma dispersion function is defined as 

exp(-t 2 ) 
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The function Z{Q) and the product £ • Z{Q) are plotted in Fig. [2] In the limiting case of M = 0, 
that is (vd — Uj) — > 0, Eq. ^ reduces to static electron and ion shielding, which corresponds to 
the spherically symmetric Yukawa potential Eq. (jl]). For |vd — Uj| ^> vt { the product Q ■ Z((^i) 
approaches —1 such that the ion screening contribution vanishes and only electron screening 



remains in Eq. see right panel of Fig. 
convergence of the Fourier integral Eq. ( 5 ) 



The finite electron contribution ensures the numerical 




2.2. Numerical implementation 

The accuracy of the Dynamical Screening Model depends crucially on the 3D Discrete Fourier 
Transform (DFT), Eq. A sufficient sampling rate in k- and r-space at the same time, requires 
a large number of sample points, which are, however, in 3D typically limited by computer 
resources to N x = N y = N z = 256. Even with the optimal choice of the sampling interval, 
corresponding to the grid spacing in the reciprocal space, the result for many of the relevant 
plasma parameters is not satisfying and lacks accuracy. In particular, the cut off error in the 
streaming direction leads to pseudoperiodicity effects from the DFT, such as artificial oscillations 
in upstream direction. 

A straightforward approach to reduce the numerical complexity of the Fourier integral 
Eq. ([5]) is to make use of the axial symmetry in streaming direction by introducing cylindrical 
coordinates. However, the numerical evaluation of the resulting radial Hankel transform is very 
time-consuming and converges poorly due to an oscillating Bessel function in the integrand. 
Instead, an alternative approach was used here to improve the memory and time efficiency by 
exploiting the radial symmetry. The 3D integration area was sliced along /^-direction and N z 
2D-DFTs were performed on a fe Xj y = 256 x 256 grid. Notably, for each /c 2 -value only the 
resulting N x /2 = 128 values in one radial direction (we used the positive half of the x-axis) 
have to be stored and used for the subsequent /^-integration. Hence, instead of N x ■ N y = 256 2 
only N x /2 = 128 DFTs in fc z -direction have to be computed. Moreover, using this scheme 
the integration area can be easily adjusted to asymmetric spatial dimensions. In particular, 
for M > 0.5 in the case of weak (Landau and/or collisional) damping, we used N z = 512 
instead of N z = 256 grid points which considerably reduced pseudoperiodicity effects and ensured 
convergence of <3?(r). 

For the sake of completeness, we should mention that because of different length scales of the 
Debye potential and the wake modulations we did not transform ^(k) directly, but the difference 
A<£(k) = c &(k) — <I>Yuk(k) [46j, where $Yuk(k) denotes the collisionless static case. Subsequently, 
utilizing the linearity of the Fourier transform, the Yukawa potential Eq. <^ is added in real 
space again. In real space, the DFT yields a grid resolution of A/) e /8 for the wake contribution 
A<£(r). The data is post processed with a spline interpolation. 

In order to ensure the convergence of the Fourier integration, sufficient collisional and/or 
Landau damping is required. More specifically, for the collisionless case (z^ n = 0) the temperature 
ratio is limited to T e /Tj < 25. For reliable results over the full range of M and temperature ratios 
up to T e /Ti = 100, weak collisional damping z^ n = 0.1uj pi is required. The Fourier integration 
performs better for small values of M. 



2.3. Results 

The LR calculations discussed here use the plasma parameters (electron temperature, density, 
ion mass, etc.) as introduced in section Neglecting the dust drift, i.e. |v^| = 0, 

the only dust parameter which enters the calculations is the grain charge, which is set to 
Qd = — 10 4 eo = —1.602 • 10~ 15 C. Since the wake potential of a point-like grain (a/A£> e — 0) 
scales linearly with Qj, a pre-computed potential <3?(r) can easily be adapted to any other grain 
charge required, see Eq. 

Fig. [3] shows the shape of the dynamically screened (wake) potential $(r) in real space 
for a range of typical plasma parameters taken from the experiments [8] 123] . To capture the 
general trends, we consider three characteristic electron-to-ion temperature ratios T r = T e /Ti = 
10,30,100 (top to bottom row) and drift velocities in the subsonic regime (M = 0.5, left), at 
Bohm speed (M = 1.0, middle column), and supersonic ion flow (M = 1.5, right). The ion- 
neutral scattering frequency is fixed at a value of Vi n /uj pi = 0.1. According to Eq. ([I]) this 
corresponds to pressures of p = 31 Pa, 18 Pa and 10 Pa for T T = 10,30,100, respectively. The 
plasma flow gives rise to an oscillating wake structure downstream the grain. On that account, 
the grain potential $(r) essentially deviates from the static Yukawa potential in all considered 
cases. Typically, the potential's anisotropy, i.e. the symmetry breaking, increases with M. The 
wakefield is most pronounced for large values of T r . 

For M = 0.5 and T e /Tj = 10 (p = 31 Pa), the range and height of the wakefield oscillations are 
strongly reduced by the overlapping effect of Landau and collisional damping. Only a relatively 
weakly pronounced primary potential maximum, i.e. a single node directly behind the grain, 
emerges in the subsonic regime j^] For M = 1.0, and even more pronounced for M = 1.5, the 
wave front becomes cone-shaped. In the supersonic case, the peak is lower than for M = 1.0. The 
same is observed for the second (third) potential peak for T r = 30 (T r = 100). Further maxima 
and minima besides the first downstream potential maximum start to appear for temperature 
ratios T r > 10. Considering T r = 100, there are in total three (significant) positive space charge 
regions, which may result in an attractive (non-reciprocal) force between the equally charged 
dust grains. For M = 1.0 the form of the wake extends far in the cross-streaming direction. 
With growing M the wake becomes increasingly stretched in the streaming direction while the 
peak height goes down. A plane $ = wave front is observed for M ~ 0.875 (T r = 30). 

The potential profile can be studied in more detail in Fig. |4j where cuts through the potential 
surface are plotted for M = 0.375,0.6875,1.0 and T v = 30. For the considered grain charge 
°f Qd = — 10 4 eo the potential height takes a value around 10~ 2 • k^T e /eo ~ 25mV, which is 
consistent with Refs. |22| 142] , Evidently, an oscillatory wake is not only formed when the speed 
of the ion flow exceeds the ion-acoustic velocity |54| 171] . Rather, a significant ion focus with 
a potential peak height of 21.7mV is already found for M = 0.375, i.e., well below ion sound 
speed. In fact, the highest amplitude of the wakefield is reached at M = 0.6875, see also Fig. [5] 
In the upstream direction, the range of $ par a is already significantly enhanced for M = 0.375, 
while in the cross-stream direction $ or tho is still relatively close to the corresponding Yukawa 
potential for that value. When M is increased, the (predominantly repulsive) potential profile 
in the upstream and cross-stream direction approach each other, as discussed in the context of 
Eq. (§. 

In the following, we will discuss the functional characteristics of the wakefield downstream 
from the grain. The potential peak height and its position as a function of M are plotted 
in Fig. [5] The positions of the maxima and minima of the wake potential are found to shift 
linearly away from the grain, when M is increased. Averaging over the three temperature ratios 

3 Compared to a isotropically screened Yukawa system, even a relatively shallow wake, such as observed for 
M — 0.5 and T r = 10, was recently shown to have a drastic effect on the collective many particle behavior [52] . 
Note that for a subsystem dust grains, the non-reciprocality along the streaming direction causes a violation of 
Newton's law of that action equals reaction, see e.g. |48| . 
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Figure 3. Linear response results: contour plots of the grain potential $(r) for varied 
streaming velocities M = 0.5,1.0,1.5 (from left to right) and electron-ion temperature ratios 
T e /Ti = 10,30, 100, that is T { = 3000 K, 1000 K, 300 K (from top to bottom). The dust grain is 
placed at the origin and the plasma flows in the positive z direction. The ion-neutral collision 
frequency is fixed to Ui n /aj pi = 0.1. Dashed-blue (solid-red) contours denote areas of negative 
(positive) space charge; potential values larger than 20 mV are clipped and shown black. The 
dotted line (on green ground) corresponds to <I> = 0. The equipotential lines are separated by 
2mV. 



T e /Tj = 10,30,50 the best linear fit a^M + bi for the i-th peak is obtained for the parameters 
Oj = (0.848,3.91,7.19, 11.1, 15.8) mm and bi = (0.0553, -0.327, -0.745, -1.54, -2.78) mm, see 
dashed gray lines in Fig. [5] (top). Corresponding PIC results for the wake maxima (see Table [T] 
in section [3| are marked with diamond symbols. Despite the finite, and relatively large grain 
size considered in the PIC simulations, there is a good agreement of both methods. On the basis 
of that result, a longitudinal wavelength of the wake equal to 27rAz) e M, discussed in Ref. |4"2] . 
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Figure 4. LR results for the electrostatic grain potential in the streaming direction z at 
x = y = (red) and perpendicular to the flow at z = (blue) for M = 0.375,0.6875, and 1.0 
(left to right panel). The dashed black lines indicate the (isotropic) Yukawa potential, Eq. Q, 
for the corresponding total Debye length Xd = 0.8452 mm. The light gray curves represent slices 
through the potential surface in z-direction for y-values in the range 2Xd s < y < 0. Considered 
parameters: T r = 30 and v% n = 0.lui Pi (cf. middle row in Fig. [3]). Note the different scaling of 
the abscissas; the range of the ordinate is [— 35 mV : 35 mV]. 

holds in the relevant range of M as a reasonable approximation. 

In Fig. [5] we also show the magnitude of the potential peak heights as a function of M. As it is 
discussed in the context of Fig. |4j the wake amplitude reveals a clear maximum when M is varied. 
Interestingly, this maximum shifts to the subsonic regime when the ratio T e /Ti is increased (in 
Fig. [5j an arrow indicates the spline interpolated maximum). According to the OML theory |72] . 
the dust charge Qqml i ncreases monotonically with M as shown for three temperature ratios 
in the bottom left panel in Fig. [5] (see the dashed black curves). Taking the charge increase into 
account, the maximum value of the peak height shifts to slightly higher values of M as displayed 
for the first potential peak (black solid line). 

The left panel of Fig. [6] shows the temperature dependence in more detail. For high ion 
temperatures, i.e. T r < 5, the largest peak amplitude is found in the supersonic regime, i.e. 
M > 1.5. Increase of T r leads to a monotonic rise of the potential height (the effect of Landau 
damping is reduced). The greatest slope of the peak height with respect to T r is observed for 
M = 0.5. So that in the limit of a large temperature ratio (a low ion temperature), the largest 
wake amplitude is found in the subsonic regime, that is M = 0.5, whereas for the largest Mach 
number, M = 1.5, the amplitude of the wake oscillation is considerably smaller. The maximum 
value of the first peak's height as function of M and its position are denoted by open black 
circles for T r = 10, 30, 50, 100. The spline through these data points reveals that the strongest 
first peak in the wake (the primary potential maximum) over the full range of Mach numbers 
shifts with increase of T r closer to the dust grain (top left panel). We note that a similar trend, 
as observed for the fixed grain charge (black dashed line) , is found when one includes the effect of 
the charge increase according to OML theory (not shown). The right panel of Fig. [6] displays that 
an increase of ion-neutral scattering V{ n leads to an increased damping of the wake oscillation 
with the peak positions being changed only slightly. 

So far, we mainly focused on the potential's wake. The form of the potential upstream and 
laterally from the grain can be captured by an effective screening length [331 [361 EQl [731 [71] . Fig. [7] 
displays the Debye length as function of M for different temperature ratios T r . For low ion drift, 
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Figure 5. Peak positions (top) and peak amplitude (bottom) of the wake potential for 
T e /Ti = 10,30,50 (left to right panel) as a function of the Mach number M. Red (blue) 
curves correspond to a positive (negative) space charge. The data points are interpolated by 
splines. In the top (bottom) panel, the peaks in the wake are counted from the bottom (top) 
curve up (down). The PIC results from Tab. [I] for the wake maxima are denoted by diamond 
symbols. The dashed black curves in the bottom left panel indicate the charge variation with M 
for T r = 10,30,50 according to OML theory (in arb. units.). The black solid curves (bottom) 
indicate the height of the first peak in the wake when the effect of charge variation is included. 

Ui "C v^, both electron and ions provide shielding, but ion screening typically dominates since 
T{ < T e . For zero flow, M = 0, the isotropic Yukawa potential, Eq. Q, applies with the total 
Debye length being A^> = 0.30,0.18,0.14 in units of \£, e for T r = 10,30,50, respectively. In 
this static case, the negatively charged dust grain Coulomb potential is predominantly screened 
by a positively charged ion cloud which isotropically forms around the grain. In contrast, in 
the supersonic limit M 3> 1, (spherically symmetric) Debye screening is entirely contributed by 
the electrons^] For finite M, the characteristic screening length is found to be different in the 
upstream and cross-stream direction, with the screening in the upstream direction being larger, 
see Fig. [7] It needs to be mentioned, that the characteristic screening lengths were obtained by 
fitting the analytic Yukawa potential to potential cuts through the dust grain such as shown in 
Fig. [4j However, these potential cuts do not exactly match the functional form of the Yukawa 
potential. For instance, for M = 0.375 (left panel in Fig. [4| it is found that the orthogonal slice 
slightly overshoots $ = 0, i.e. the potential possesses a weakly attractive branch sideways from 
the grainj^] Therefore, the screening lengths in Fig. [7] cover the form of the potential close to 

4 See also the discussion of Eq. jfij) in section 2.1 

The attractive part is found in a small parameter range only. In experiments and PIC simulations this weakly 
pronounced effect has not been observed so far |73| . 
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Figure 6. Positions (top) and heights (bottom) of the first peak of the wake potential 
for different M values as a function of temperature ratio T e /Tj (left panel, V{ n = 0.1) and 
collision frequency v^ n (right panel, T e /Ti = 100). The data points are denoted by symbols and 
interpolated by splines. Open black circles mark the absolute maximum of the peak height as a 
function of M (bottom left), cf. Fig. p] and its position (top left) for different ratios of T r . 



the grain, but have an approximate character only. Also, a systematic error is observed in the 
limit M — > oo, where the value of the effective screening length in upstream direction converges 
towards a value that is about five percent larger than the appropriate screening length that is the 
electron Debye length Xo e - In general, we find that the Yukawa potential matches the potential 
profile perpendicular to the flow essentially better than in the upstream direction which is in 
accordance with [73J. 

Finally, we like to point out that the characteristic screening length can be roughly 
approximated by the interpolation formula 

kBTi+rjm^ \ 1/2 / 1 , f s (l) 



Xs Xde \k B T e + k B T t + m iU j ) [mx 2 Di + mxij ' (9) 

where f s (rj) = 1 + r/M 2 ^ and rj is a heuristic fitting parameter which we introduce here. 
Considering r] = 1, A s reproduces both screening limits correctly: (i) X s = for M = 0, 
and (ii) A s = Xjj e for M — > oo (see the dotted lines in Fig. [7J [36] , However, X s is found to 
converge much slower towards Xd £ than the curves obtained by LR theory (Fig. [7J and PIC 



simulations (as will be shown in section 3.3, Fig. 11). Therefore, it is suitable to introduce the 



fitting parameter 77 = 1.27 which allows for an essentially better agreement with the present 
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Figure 7. Linear response (LR) results for the effective screening length in the upstream (solid) 
and cross-stream (dashed) direction relative to the grain as a function of the ion flow velocity M. 
The LR results are given for two temperature ratios T r = 10 (red) and T r = 50 (blue). The light 
gray [black] lines indicate the characteristic screening length \ s , Eq. ([9]), for rj = 1.0,1.27, and 
A fc , Eq. (FLO}, for T r = 10 [T r = 50], respectively. 



results in the relevant range < M < 2j^] The concept of an effective screening length in a 
plasma with streaming ions was also discussed by S.A. Khrapak et at, see |74j and on page 45 
in Ref. [5]. Here, the analytical form of the effective screening length is approximated by the 
formula 



-1/2 

\ k =[f m (M T )\ D ] + \- D 2 e ) (10) 

where the two fitting functions /i = exp(— My/2) and /2 = (1 + My) -1 are defined in terms 
of the thermal Mach number My = Uj/fy = Mc s /fy. Evidently, with the fitting function /2, 
Eq. ( 10 ) is identical to Eq. Q with r\ = 1 and is therefore not further considered in the following. 



For T r = 10, Eq. (10) yields in conjunction with f\ a very reasonable approximation (especially 
for the screening in the upstream direction) without any heuristic parameter, see dashed gray 
line in Fig. [7] Furthermore, A& shows (in contrast to A s ) the correct analytic trend of a rapidly 
increasing slope with increasing temperature ratio T r , which is evident in the intersection of 
curves at different temperatures. However, overestimates the decay of the ion contribution to 



the screening as shown for T r = 50 in Fig. [7] and for T r = 30 in Fig. 11 



Theoretical findings that the effective screening length converges toward the electron Debye 
length \£) e as M increases is supported by experimental results |75[ 176] . In those experiments 



6 At larger values of M screening by the electrons plays the dominant role and the effective screening length is 
therefore close to A_o e . 



the effective screening length for the grains suspended in the sheath was found by analyzing 
collisions between the dust grains. In both cases it was shown that supersonic ions are unable 
to screen the dust particles and therefore the electron Debye length is the appropriate screening 
length for dust grains in supersonic plasma flows. 

3. 3D Particle-In-Cell Simulations 

The LR theory, employed in the previous section, neglects some important processes associated 
with the dust grain charging and wake formation. In particular, in the regime of strong plasma 
flows and high grain charges, non-linear effects can be significant. Particle kinetic simulations, 
such as the PIC method, allow for self-consistent studies of the grain charging, and account also 
for non- linear phenomena. 

3.1. Theoretical Approach 

Analytical plasma models that account also for non-linear effects are difficult to develop. 
Therefore, one often employs numerical kinetic particle simulations, in which plasma particle 
trajectories are followed in self-consistent force fields without many approximations, allowing for 
a self-consistent evolution of the system. As plasma particle trajectories are the characteristics 
of the Vlasov equation, the PIC simulations can be considered as an alternative way of solving 
the Vlasov equation for arbitrary particle distributions |77j . 

Plasma simulations, in which forces are evaluated between all pairs of simulation particles, are 
typically very expensive in terms of computational time and memory. A simulation of n particles, 
with such direct force calculations, has complexity 0(n 2 ), as approximately n 2 equations need 
to be solved to find the forces on all the particles. O relates to the computer resources usage 
time. As the number of simulated particles needs to be large, especially in 3D simulations, such 
an algorithm would be very inefficient. 

Introducing a grid within the PIC scheme can significantly reduce the complexity of the 
algorithm |78| . In the PIC method, the physical quantities of each simulation particle are 
weighted to neighboring grid points to build appropriate density fields on the grid. In the 
electrostatic approximation, the electric potential is found from the charge density on the 
grid points by solving the Poisson equation. Thus, the complexity of the PIC algorithm is 
O(n) + 0{n g log(n g )), where n g denotes the number of grid points with n g <C n. 0{n g log(n s )) 
is the complexity of solving the field equations. 

Integration of particle trajectories puts limits on PIC simulations. With both electrons and 
ions simulated, the time resolution should be a fraction of the electron plasma period. The 
charging time is usually of the order of the ion plasma period, and thus it requires a large number 
of time steps. Simulations can be accelerated by assuming Boltzmann distributed electrons, which 
leads to a hybrid PIC-fluid code |79j . This approximation is however not always valid, as trapped 
electrons, or electron sources (e.g., due to photoemission) can give rise to local deviations from 
the Boltzmann distribution |80) . Another way of speeding up the evolution of the system is by 
using a reduced ion mass mj , as the ion plasma frequency increases with decreasing m; . Ion mass 
values as low as mi = 30m e have been used in the literature |81j . In |37| it has been demonstrated 
that mass ratios rrii/m e > 100 give reliable results for the dust charging. While the reduction 
of ion mass leads to some quantitative differences, the results are qualitatively correct, scalable 
and can be presented in normalized units. Note that with the ion mass being reduced, the sound 
speed c s is increased and the floating potential is slightly reduced. In turn, with a relatively 
large grain, the discreteness noise and charge fluctuations on the grain are diminished. 

3.2. Numerical implementation 

In the present study, we use the Dust in Plasma 3D (DiP3D) particle-in-cell code, which 
has been developed to study the charging of dust grains and other objects in various 



plasma environments |41| 15 Ij . The code simulates a collisionless plasma in the electrostatic 
approximation, with electrons and ions represented as individual plasma particles. To weight 
particle charges to the regular grid, and to project forces from the grid to the particles, first 
order linear weighting is used. The particle trajectories are advanced with the leap frog method 
characterized by a staggered time mesh [78j: 

Xi(t + At) = Xi(t) +Vi(t + At/ 2) At , 
Vi(t + At/2) =Vi(t- At/2) + f i (t)At/m l , (11) 

where i refers to an electron or ion, fi is the force projected on the i-th particle from the nearest 
grid points, and At is the computational time step. 

The plasma particles are initially uniformly distributed with Maxwellian velocity distributions 
in the simulation box of size L 3 = (24A£) e ) 3 . The particles are injected into the simulation box at 
each time step according to prescribed particle fluxes. To simulate an open plasma system, the 
particles can leave freely through the boundaries of the simulation box. The code uses Dirichlet 
boundary conditions for the potential. The potential at the boundaries is set according to the 
plasma potential = V. 

One or several spherical dust grains are placed inside the simulation domain far away from 
the boundaries. We assume a conducting grain, and each plasma particle that hits the surface 
is removed and contributes to the current to the grain. The charge of the removed particle is 
distributed on the grain surface to cancel internal electric fields. The grain is initially uncharged 
and becomes charged self-consistently during the simulation. 

The DiP3D code is run on a computer cluster, with typically n = 10 7 or more simulation 
particles distributed on several nodes, with a typical grid size of n g > 128 3 . The Message- 
Passing-Interface (MPI) is used for parallelization. The code can be stopped and then restarted 
with modified dust grain configurations, which can reduce the computation time in case of 
systematic parametric studies. The code also allows for a movement of the simulated object, 
and for the inclusion of various charging processes such as photoemission. Recently, the code has 
been upgraded to include external magnetic fields and ion-neutral and electron-neutral collisions 
|82j . With all the features, the code can be placed among other cutting-edge particle plasma 
codes for object-plasma interactions |36| 156) I83j . In the present work, we make only use of the 
basic features of the code. 

It should be noted that the dimensionality of the code is an important aspect of any simulation. 
For the ultimate goal of direct application of the numerical results to support laboratory studies, 
3D simulations should be used. Many important effects, such as ion focusing, are still present 
when the dimensionality is reduced. However, in a 2D system, the charge of the plasma particle 
corresponds to the charge of an infinite rod, and a simulated circular grain represents an infinite 
cylinder. The plasma frequency cj p , Debye length Xp, and charge density do not change with 
the dimensionality of the system. Reducing the dimensionality can give faster algorithms, but 
while the results of such simulations are qualitatively correct, they do not need to represent 
a 3D system quantitatively. Thus, results from codes with reduced dimensionality should be 
compared with corresponding theoretical models. In the present paper, we present results from 
3D simulations. 

3.3. Results 

With the DiP3D code, we performed simulations of a single grain in a collisionless plasma. The 
grain radius is a = 0.185A/) e = 0.1563 mm. The finite size of the grain leads to enhanced potential 
variations in the wake as compared to the previous LR results, because the self-consistent grain 
charge is on the order of Qd ~ — 10 5 eo- All other plasma parameters are as described in 
section [T7T] The grain charging and wake formation have been investigated for different electron- 
to-ion temperature ratios, T r £ {10, 30, 100}, and Mach numbers, M £ {0.75, 1, 1.5}, as shown in 
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Figure 8. First principle results from collisionless 3D PIC simulations. The contour plots 
show the wake-potential <3? in the y-z plane at x = for varied ion streaming velocities 
M = 0.75,1.0,1.5 (from left to right) along z axis, and electron-ion temperature ratios 
T r = 10,30,100 (from top to bottom). The finite-sized dust grain is placed at the origin. Blue 
(red) color codes denote areas of negative (positive) space charge; potential values larger than 
= 0.2 V (lower than $ = —0.2 V) are clipped and shown black (white). The size of the dust 
grain is marked by the black, open circle at the origin of the coordinate system. 

the contour plots of the electric potential in Fig. [8] In all cases, we observe a positive maximum in 
the potential distribution in the wake at a distance of approximately one to two \d £ downstream 
from the dust grain. This potential maximum is related to the ion focus in the wake of the 
negatively charged grain. The maximum in the potential is located further downstream from 
the grain than the enhancement in the ion density in the focus region, which is consistent with 
the Poisson equation V 2 $ = — p/eo, where p is the charge density. Ion focusing is due to the 
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rel. position \Xo e = 845.2^/n] rel. position [A/^, = 845.2jlm] rel. position \Xo e = 845.2/iml 

Figure 9. 3D PIC results for T e /Tj = 30: potential along the flow direction z through the 
centre of the grain x = y = (red line) and in the cross-stream direction at z = (blue line) 
for M =0.75, 1.0, and 1.25. Light gray curves represent slices through the potential surface for 
< y < 2A£> e (along the direction of the flow). The dashed black line shows the potential drop 
for the static case, M = 0, for which = —2.02V. As the floating potential of the grain changes 
with the flow speed, the results for the static case have been rescaled, such that the potential on 
the grain matches the floating potential on the grain in a flowing plasma. 



deflection of ion trajectories by electric fields in the vicinity of a negatively charged grain 

Similarly to the LR results in section [2] in the PIC simulations the potential enhancements in 
the wake become more pronounced for larger M and larger electron-to-ion temperature ratios T r , 
i.e., colder ions. The first local potential maximum is followed by other potential extrema further 
downstream. For T r = 100 and M = 0.75, the pattern in the wake shows a diverse structure with 
several minima and maxima being located also away from the flow symmetry axis. Note that this 
could also be an artefact due to the Poisson solver on the finite grid. In particular for very low M 
and T r , the possibility of instabilities can become an important issue for the convergence of the 
code [84j. For larger M and T r , the extrema in the potential distribution are along the direction 
of the flow. With an increase of the ion temperature, the Landau damping becomes important, 
and thus ion focusing and wake patterns are less pronounced. For T r = 10 and M = 0.75 only 
a weak, single maximum is observed. The considerable structural consequences of even a weak 
ion focus were subject of recent investigations |52j . 

In Fig. [9j we present potential cuts along the flow direction z at y = and x 6 [0, 2Au e ] 
for T r = 30. In the figure, we also plot the potential cut in the cross-stream direction, as 
well as the potential cut for the grain in stationary plasma. In contrast to the LR results, the 
potential maximum monotonically increases with the flow speed, and so does its extension in the 
perpendicular direction whilst the Mach cone forms. However, there is an important difference 
between the LR and PIC results. In the LR analysis, the charge on the grain is constant, while 
in the PIC simulations it changes with the flow. The floating potential on a dust grain in a 
stationary plasma is $g = —2.02 V and it becomes larger with increasing M (see discussion 
below). While the charge and potential on the grain increase with the flow, the amplitude of 
peaks in the wake potential is being enhanced. 

In order to compare the PIC results with LR calculations, in the left panel of Fig. 10, we 
present potential cuts through the grain along the z-direction and in the cross-stream direction 
for M = 1.0 and T e /Tj = 30 as obtained from collisionless 3D PIC simulations and from LR 
calculations with collisional damping included (i / i n /u) Pi = 0.1). The Yukawa potential, Eq. 
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Figure 10. Left: Potential cuts through the grain along the z-direction (solid lines) and 
in cross-stream direction (dash-dotted lines) for M = 1.0 and T e /Tj = 30 as obtained from 
collisionless 3D PIC (red) and LR calculations with Vi n /oj pi =0.1 (blue). The grain charge of 
the LR result and the corresponding Yukawa potential are adjusted to match the potential of the 
PIC simulation, which includes the grain charging process self-consistently. Note that in contrast 
to PIC, the LR results do not take into account the finite particle radius a = 0.185A£> e , which 
explains the offset of the curves. Right: Potential cuts along the z-direction for M = 1.0 and 
T e /Ti = 50. Shown are our LR calculations with finite damping included Vi n /uj pi = 0.1 (blue 
solid line), COPTIC data by Hutchinson for a point-like grain in a collisionless plasma (black 
solid line) |42], and linear response results for the collisionless case by Lampe et al. (red dotted 
line) |33] , 

for a grain in stationary plasma is also shown. The grain charge of the LR result is adjusted 
to match the potential of the first peak in the PIC simulation (the same charge is used for the 
Yukawa potential). It can be inferred from the plot that the agreement between the PIC and 
LR results is good for the first two extrema. It is important to note that in contrast to the PIC 
simulations, the LR results do not take into account the finite particle radius a = 0.185A£> e - For 
this reason, the LR curves are slightly shifted towards zero. When the finite size of the grain is 
taken into account, the agreement between the LR and PIC results for the potential profiles is 
improved. 

As discussed in the context of Fig. [5l we find an overall good agreement in the longitudinal 
wavelength of the PIC and LR results^] A more inconsistent discrepancy in the wavelengths 
was recently reported by Hutchinson [42 J. Therefore, in the right panel of Fig. [l0]a comparison 
for T e /Ti = 50 and M = 1.0 is shown between our LR results, COPTIC results for a point-like 
particle (Fig. 11 in [42J), and the linear calculation by Lampe et al. (Fig. 3 in |33|). The main 
difference between Lampe et a/.'s and our LR results is due to the fact that we include finite 
collisional damping (which in our case is necessary to avoid pseudoperiodicity effects for the given 

7 In accordance with |35j . the linear approach tends to overestimate the asymmetry of the potential distribution 
around the dust particle. 
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Figure 11. Effective screening length of the dust grain potential for T r = 30 as a function 
of the ion flow velocity M. PIC [LR] results for the upstream (solid curve) and cross-stream 
(dashed) direction are indicated by circles [crosses] on orange/red [black] lines. These lines are 
result of an approximating spline interpolation and therefore do not exactly pass through the 
data points. The estimated error of the PIC data points is 5 < 0.05\D e - The light gray lines 
indicate the screening length A s , Eq. (|9]), for rj = 1.0,1.27, and Xk, Eq. (10), respectively. See 
also Fig. [7] 



temperature ratio T e /Tj = 50, as stated in section 2.2 ). The wavelength obtained from COPTIC 
is in full accordance with our LR calculation]^] The stronger decay of the wake oscillations can be 
attributed to finite collisional damping which is included in the LR calculations (vi n /bJ Pi = 0.1). 
However, in the left panel of Fig. [5] we find a much better agreement in the wake amplitude 
at the first subsidiary valley. The stronger damping in our PIC simulation (left panel) must be 
due to the large grain size (and charge) giving rise to much stronger non-linear effects than the 
point-like grain in the COPTIC simulation. As discussed in |l2], a major effect of non-linearity 
is the local increase of the effective ion temperature in conjunction with an enhanced Landau 
damping of the wake amplitude. The presence of (strong) non-linearities in the present PIC 
(DiP3D) simulations becomes apparent by the similarity of the shape of the primary potential 
maximum for M = 1 and the temperature ratios T r = 30 and especially T r = 100 (Fig. [8]) and 
the non-linear wake structure presented in Fig. 5(a) in Ref. |42] , Therefore, we assume that the 
good agreement between our collisionless PIC simulation and the LR calculation with collisions 
included — not only in the wavelength of the wake oscillations but also in the decay of the wake 
amplitude — is due to a nonlinearity-induced damping mechanism that reduces the amplitude of 
the wake potential. 

In the downstream direction, the primary repulsive branch of the potential is shifted closer 



Collisional damping has on a minor effect on the peak positions, see right panel in Fig. 6 



towards the grain at finite (subsonic) flows, which is similar to the LR results (cf. Figs. [4] and [9]). 
As ions contribute more to the grain screening in stationary than in flowing plasma, the grain 
screening is strongest in the static (Yukawa) case and weakens for finite M, see again Fig. [9] The 
quantitative comparison of the PIC and LR results has been done by calculating the effective 
screening length from the potential profiles for different M for T r = 30 in the cross-stream and 



upstream directions, see Fig. 11 The general trend for the PIC and LR results is that the 
effective screening length A increases with the flow velocity M. The effective screening for the 
stationary plasma in the PIC simulations is A = 0.23A£> e . This value is larger than the the results 
from the LR calculations, as well as larger than the total Debye length, which is Xrj = 0.18A£> e . 
The difference is due to a finite grid resolution in the simulation, where the grid spacing is 
Ax ~ Xd- For the flow velocities M < 0.5, the cross-stream and upstream shielding lengths 
are lower than 0.4A/) e . The values are roughly 0.1\r) e larger (smaller) than the cross-stream 
(upstream) results from LR calculations. For M > 0.5, the effective screening length from the 
PIC simulations strongly increases with M, with steeper increase for the upstream direction. 
This is the same behavior as for the LR results. At M < 1, for both LR and PIC results, the 
screening length A is larger on the upstream than on the downstream side. The increase of the 
effective screening length with the flow is in agreement with previous studies |33| 170] . however, 
the details of the trend for cross-stream and upstream directions are different. In Ref. |33j the 
cross-stream shielding was systematically larger than the upstream one. Based on our results 
from simulations with two different numerical methods, the upstream shielding length is larger 
than the cross-stream one for M < 1. For higher M, the trend continues for the LR calculations, 
while for PIC simulations, the shielding becomes larger on the cross-stream than on the upstream 
side. The effective screening length converges towards the electron Debye length for supersonic 
flows, as it was discussed in the previous section for the LR results. 

Finally, in Table [T] we summarize the floating potential on the grain as well as the positions 
and heights of the maxima for different M values. Since the ion flux to the grain surface changes 
with the ion flow speed, the floating potential <3?y/ on the grain becomes more negative with 
increasing M. Changing T r has little influence on the floating potential, as <3?fl is mainly controlled 
by the electron temperature T e , which is kept constant in our simulations. The heights of the 
potential maxima increase with the flow speed, and so does the distance between the peak and 
center of the grain. 



4. Method-Spanning Comparison 

We now summarize the correspondence of the two methods and discuss their limitations. 



4-1- Comparison of the potentials 

The dynamical screening potentials obtained from the LR approach are topologically similar to 
the corresponding results from the PIC simulations, see again Figs. [3] and [8] The PIC results 
are more noisy, as they account for the dynamics of plasma particles and fluctuations in the 
system. Since the PIC simulations consider a finite-sized grain that is self-consistently charged, 
the potential on the grain changes with the flow speed, see Table [T] The comparison of the peak 
positions of the positive ion wake charges, Fig. [5] shows a good agreement. Note that the exact 
position of weak peaks with a low signal to noise ratio (high peak numbers at small T r and M 
values) can just barely be resolved, cf. Fig. [3] 

The amplitude of the peak in the wake calculated with LR achieves a maximum at a certain 
ion flow velocity (Fig. [5]). This maximum is at lower flow velocities for larger temperature ratios. 
We do not see a maximum in the peak amplitude for the PIC results summarized in Table [T] 
However, in section [2j we considered the charge on the grain to be constant, while in the PIC 
simulations the potential on the grain changes with the flow. When scaling linearly all PIC 
results to a fixed potential, we obtain a maximum in the height of the peak as a function of flow 



Table 1. Floating potential on the grain and the position pj and height hi of the i-th 
maximum in the wake. For some lower velocities only the first maxima were observed. Crosses 
(x) correspond to possible maxima located outside the simulation domain. The floating potential 
on a grain in a stationary plasma is <3?fl = —2.02 V (T r = 30). Note that these results apply to 
the reduced ion mass. (In particular, the floating potential is therefore reduced, as discussed in 
section 3.1 ). 
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velocity. This maximum is at lower velocities for higher temperature ratios, similarly to the LR 
results. Moreover, adjusting the correct grain charge, the potential pattern in the wake from the 
LR simulations agrees well with the PIC results, see again Fig. [10| 

Both, LR theory and PIC simulations confirm that the effective screening of the grain changes 
with the flow. While the effective screening length calculations with the PIC method can 
be considered as less accurate due to potential representation on a finite grid, the trend for 
the effective screening length agrees with the LR results. For subsonic ion flow velocities, the 
screening length on the upstream side is usually larger than on the cross-stream side. 

4-2. Consideration of shadowing effects 

An important question for the OCP model is to determine to what extent the wake potential 
around two or more grains can be approximated by the linear superposition of single grain 
potentials. We will thus study the impact of non-linear effects in the multiple grain wakes by 
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Figure 12. (a) Potential around two grains separated by approximately Xd b , for T e /Tj = 30 and 
M = 1.25 obtained by PIC simulations, (b) The potential constructed by a linear combination of 
the single grain wake potentials, (c) The potential difference between the constructed potential 
and the results from the simulation = $ii near — <£. Note the different values on the color 
bars. 
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Figure 13. Same as Fig. 12, but for two flow-aligned grains (M = 1.25). 



comparing the results from PIC simulations of two grains with the results obtained by a linear 
combination of data from the corresponding simulation of a single grain. 

In Fig. (12^ a), we show results from the PIC simulations of two grains for T e /Xi = 30 and 
M = 1.25. The downstream grain is located off the symmetry axis at a distance 0.93Xr> e 
from the upstream grain. The charge of a downstream grain is only affected a little by the 
upstream grain, as the downstream grain is located out of the ion focus region |51j . In the linear 
combination of two wakes originating from a single grain calculation, the grain charges (and 
thus the grain potentials) are adjusted in such a way that they match the corresponding values 
from the simulations of two grains, see Fig. 12 'b). The constructed wake captures the general 
features of the results from the simulation of two grains, i.e., the position and approximate height 
of the potential maxima and minima. However, the small differences between the two results 



indicate that non-linear wake effects are present in the simulations of two grains. In Fig. 12 




(c) the potential difference A$ = $ii nea r — & between the linear combination of wakes and the 
original result is plotted. In the system of two grains, the downstream grain influences the ion 
dynamics in the wake, and therefore leads to an asymmetric ion focus region, which can not be 
fully represented by a linear combination. The potential resulting from the linear combination 
overestimates the maximum in the wake by up to 20%. Also the symmetry of the wake potential 
is broken. 

The corresponding results for two grains aligned with the flow that are separated by 0.74X£> e 
are shown in Fig. 13 The charge on the downstream grain is strongly affected by ion focusing, and 
it is only about 30% of the charge on the upstream grain. Thus in the linear combination of the 
wake, the contribution for the downstream grain has been adjusted accordingly. The constructed 
wake again represents well the topology of the wake from the simulations of two grains. In this 
case, the relative potential difference is up to 30% with the strongest discrepancies associated 
with the first maximum. 

The agreement between constructed and simulated wake is worse for subsonic flows. Figs. 



14 and 15 show corresponding results for the subsonic case M = 0.75 for the downstream grain 



located off axis and aligned with the flow, respectively. In both cases, the differences A$ are up 



to 50% in the closest vicinity of the grains. 

In all considered cases the general features of the wake, such as positions of the maxima, are 
well reproduced. In the constructed wakes, the heights of the first maximum are overestimated, 
and for asymmetric grain arrangements, the symmetry of the wake further downstream from the 
grains is also modified. This is due to a single asymmetric ion focus region forming in a system 
of several grains with small separation which are not aligned with the flow. Thus, constructing 
the wake from the linear combination of two wakes of single grains should be taken with care. To 
construct such a wake it is necessary to know a priori the charge or potential on a downstream 
grain, which can be strongly influenced by the wake effects |50[ 151] , While the main features of 
the wake can be well represented by the linear combination of a single-grain wakes, there may 
be inaccuracy due to non-linear effects associated with ion focusing. 

4-3. Range of applicability of both simulation models 

In order to study the wake structure behind a grain in streaming plasmas, we have employed two 
complementary methods: (i) the numerical evaluation of the linearized electrostatic potential 
that is created by a point-like charge, and (ii) self-consistent, full 3D particle-in-cell simulations. 
These methods have the best performance in different regimes. The LR method has the best 
convergence for small T r and large V{ n (and of minor importance small M), i.e., large Landau or 
collisional damping are optimal for the performance of the code. The PIC simulations are easiest 
for higher M and small T r , as the simulations for small M and cold ions are often favorable 
condtitions for instabilities [84] . 

In the full PIC method, the ions have often reduced mass. The results with reduced ion 
mass are reliable, when the analysis is carried out in dimensionless units |37]. The PIC method, 
which includes self-consistent charging of the grain, accounts for a finite grain size, while the LR 
assumes a point-like dust grain with a given charge. The considered particle radius a = 0.185Xo e 
in the PIC simulations corresponds to the non-linear regime [42j. 

The two approaches are complementary, and together allow for studies of a grain potential over 
a wide regime of parameters. Moreover, for single grain simulations, there is good agreement in 
the functional form of the wake (i.e., number and relative positions of the wake potential maxima 
and minima). The main issue of the LR approach is the correct estimate of the grain charge 
as a function of M, and also for a grain that is shadowed by another grain located upstream. 
A possible solution to overcome this issue could be the use of pre-computed tables obtained 
from self-consistent PIC simulations or making use of data from experiments |24j . It is also 
worth mentioning that possible deviations from a shifted Maxwellian ion velocity distribution, 
as reported in |42 ^ 85 | l86]. need to be considered in developement of advanced models of the dust 
plasma interactions applicable for a wide range of plasma parameters. 

5. Summary and Outlook 

In the present work, we have investigated the electrostatic potential distribution around a 
spherical dust grain in stationary and flowing plasmas by means of linear response calculations 
and particle-in-cell simulations. The streaming plasma leads to strong deviations from the 
statically screened Coulomb potential, giving rise to a distinct oscillatory wake structure behind 
the grain. For both sub- and supersonic flow velocities, positive ion wake charges downstream 
from the grain are present that can give rise to non-reciprocal dust grain interactions. We 
have presented a systematic overview of the LR results for different flow velocities, electron-to- 
ion temperature ratios, and ion neutral collision frequencies. The LR results have been tested 
against self-consistent particle-in-cell simulations. The results from the PIC simulations for the 
wake pattern agree qualitatively with the LR data when the finite size of the grain is accounted 
for. 



Downstream from the grain, a pronounced oscillatory wakefield with several maxima develops 
with increasing ion flow. Independent of the temperature ratio T e /Tj, the peak positions are 
shifted linearly away from the grain when M is increased. Upstream and to the side from the 
grain, the potential decay can be approximated by a Yukawa potential. The effective screening 
length increases in both directions with the flow, but in contrast to earlier LR results [33j, it is 
found to be larger in the upstream direction than laterally from the grain. 

In a system of several grains in a plasma flow, the charge and potential on the downstream 
grains can be significantly modified by wake effects. The results from PIC simulations show 
that the wake behind the two grains can include non- linear effects, and that the charge on 
the downstream grains can be significantly reduced. While the general wake pattern can be 
represented by a linear combination of wakes behind a single grain to a rather good accuracy, 
local deviations up to 50% of the potential are observed for subsonic flows in the closest vicinity 
of the grains. Thus, special care needs to be taken when applying LR results in OCP simulations, 
and at least the charge reduction on the downstream grains needs to be considered for the model 
to be credible. 

The influence of wake effects on collective many-particle behaviour in 3D plasma crystals has 
become a very timely question and is currently under ongoing experimental and theoretical 
investigation |43[ HU HSJ [52]. The presented wake potentials are currently employed in 
Dynamical Dust Simulations which accompany the experimental exploration of string formation 
and externally triggered phase transitions in 3D confined dusty plasma clouds |44t 145] , OCP 
simulations that include many particles require, however, an adequate adjustment of the charges 
of the individual grains, which are shadowed by the upstream grains. Such data may be obtained 
from the first principle PIC simulations or from experiments [24J. The development of such a 
grain-charging model is the subject of ongoing work. 

Due to the high scalability of Coulomb systems, we expect these results to be of interest also 
to other types of multi-component plasmas. In particular, the generalization of the dynamically 
screened dust approach will open the way to a description of non-equilibrium quantum systems, 
such as warm dense matter, on the microscopic level. There, the classical dielectric function 
must be replaced by the corresponding quantum (Mermin) dielectric function |54 [ I87 [ [88] . This 
scheme then allows for the computation of collective many-body properties of strongly correlated 
ions embedded into a partially ionized quantum plasma of electrons by first principle molecular 
dynamics simulations. 
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